Purpose

  • Detect start and end of periods in a shroeder using amplitude cross-correlation
  • Figure out a way to measure fine scale structural variation and periodicity in schroeders

 

   

wave_col <- viridis(10)[7]

Exploring data

sound file SchroederCallsTest4_3Nov2022.wav:

schroeders <- read_wave("./data/raw/SchroederCallsTest4_3Nov2022.wav")

spectro(schroeders, scale = FALSE, osc = TRUE, fastdisp = TRUE, grid = FALSE,
    palette = viridis, collevels = seq(-100, 0, 5), wl = 512, colwave = wave_col,
    heights = c(1, 1))

First bout

Spectro + oscillogram

bout_1 <- cutw(schroeders, from = 0, to = 0.2, output = "Wave")

spectro(bout_1, scale = FALSE, osc = TRUE, fastdisp = TRUE, grid = FALSE,
    palette = viridis, collevels = seq(-100, 0, 5), wl = 512/2, colwave = wave_col,
    heights = c(1, 1))

Oscillogram

oscillo(bout_1, colwave = wave_col)

Oscillogram zooming in

oscillo(bout_1, from = 0, to = 0.03, colwave = wave_col)

Amplitude cross-correlation

Pearson correlation

Try it on the first bout. First get correlation for an amplitude template with the first few samples

# make a template from the first fifth of duration of the entire
# shroeder
template <- bout_1@left[seq_len(length(bout_1@left)/5)]

cors <- vapply(seq_len(length(bout_1@left)/4), function(x) {

    segment <- bout_1@left[x:(x + length(template) - 1)]

    cor(template, segment)
}, FUN.VALUE = numeric(1))


plot(cors, type = "l", col = "orange", xlab = "samples", ylab = "Pearson correlation")

Then get location of correlation peaks

tpks <- fpeaks(cbind(1:length(cors), cors))

# get highest peak
tpks <- tpks[tpks[, 2] > 0.8, ]

# order by sample number
tpks <- tpks[order(tpks[, 1]), ]

# remove first peak
tpks <- tpks[-1, ]

Then get the mean difference (in samples) between peaks (which is inversely related to the frequency)

sample_diffs <- diff(tpks[, 1])

mean_diff <- mean(sample_diffs)

print(round(mean_diff))
## [1] 244

Now get the start of each period in the sound clip

# get envelope peaks
amppk <- fpeaks(cbind(1:length(bout_1@left), abs(bout_1@left)/max(abs(bout_1@left))),
    plot = FALSE)

# filter highest peaks
amppk <- amppk[amppk[, 2] > 0.95, ]

# get the highest one
first_max <- amppk[which.max(amppk[, 2]), 1]/bout_1@samp.rate

# plot
oscillo(bout_1, from = 0, to = 0.03, colwave = wave_col)

abline(v = first_max + (-100:100 * (mean(sample_diffs)/bout_1@samp.rate)),
    col = "red", lty = 2, lwd = 2)

Create a function to spot periods

spot_periods <- function(wave, plot = TRUE, abs.amp = TRUE, alpha = 0.5,
    from = 0, to = 0.03) {

    # make a template from the first fifth of duration of the
    # entire shroeder
    template <- wave@left[seq_len(length(wave@left)/5)]

    cors <- vapply(seq_len(length(wave@left)/4), function(x) {
        segment <- wave@left[x:(x + length(template) - 1)]

        cor(template, segment)
    }, FUN.VALUE = numeric(1))

    tpks <- fpeaks(cbind(1:length(cors), cors), plot = FALSE)

    # get highest peak
    tpks <- tpks[tpks[, 2] > 0.8, ]

    # order by sample number
    tpks <- tpks[order(tpks[, 1]), ]

    # remove first peak
    tpks <- tpks[-1, ]

    sample_diffs <- diff(tpks[, 1])

    # mean distance in samples between amplitude peaks
    mean_diff <- mean(sample_diffs)

    norm_abs <- if (abs.amp)
        abs(wave@left)/max(abs(wave@left)) else wave@left/max(wave@left)

    # get envelope peaks
    amppk <- fpeaks(cbind(1:length(wave@left), norm_abs), plot = FALSE)

    # filter highest peaks
    amppk <- amppk[amppk[, 2] > 0.95, ]

    # get the highest one
    first_max <- amppk[which.max(amppk[, 2]), 1]/wave@samp.rate

    # get starts of periods
    positions <- first_max + (-1000:1000 * (mean_diff/wave@samp.rate))

    # filter positions within duration of wave
    positions <- positions[positions > 0 & positions < duration(wave)]

    # plot
    if (plot) {
        osc <- oscillo(wave, from = from, to = to, colwave = wave_col)

        abline(v = positions, col = adjustcolor("blue", alpha.f = alpha),
            lty = 2, lwd = 3)

        points(y = wave@left, x = seq_len(length(wave))/wave@samp.rate,
            type = "l", col = wave_col)
    }
    return(positions)
}

Test function

Try function on other shroeders

5th bout

bout_5 <- cutw(schroeders, from = 4.8, to = 5, output = "Wave")

pos <- spot_periods(wave = bout_5, alpha = 0.5)

9th bout

bout_9 <- cutw(schroeders, from = 10.2, to = 10.4, output = "Wave")

pos <- spot_periods(wave = bout_9, to = 0.02)

11th bout

bout_11 <- cutw(schroeders, from = 15, to = 15.2, output = "Wave")

pos <- spot_periods(wave = bout_11, to = 0.02)

20th bout

bout_20 <- cutw(schroeders, from = 31.8, to = 32, output = "Wave")

pos <- spot_periods(wave = bout_20, to = 0.015)

Dynamic-time warping

# make a template from the first fifth of duration of the entire
# shroeder
template <- bout_1@left[seq_len(length(bout_1@left)/8)]

dists <- pbapply::pbsapply(cl = 4, seq_len(length(bout_1@left)/4),
    function(x) {

        segment <- bout_1@left[x:(x + length(template) - 1)]

        dtw_dist <- dtw::dtwDist(mx = rbind(template, segment))
        return(dtw_dist[1, 2])
    })

saveRDS(dists, "./data/processed/dtw_distance_bout_1.RDS")
dists <- readRDS("./data/processed/dtw_distance_bout_1.RDS")

dists <- dists/max(dists)
sims <- 1 - dists

plot(sims, type = "l", col = "orange", xlab = "samples", ylab = "DTW similarities")

Then get location of correlation peaks

tpks <- fpeaks(cbind(1:length(sims), sims), threshold = 0.3)

Making the time series stationary (diff())

stat_sims <- diff(sims)
stat_sims <- stat_sims + abs(min(stat_sims))
stat_sims <- stat_sims/max(stat_sims)

tpks <- fpeaks(cbind(1:length(stat_sims), stat_sims), threshold = 0.98)

# get highest peak
tpks <- tpks[tpks[, 2] > 0.98, ]

sels <- data.frame(sound.files = 1, selec = 1:nrow(tpks), start = tpks[,
    1], end = tpks[, 1] + 6, peak = tpks[, 2])

# merge peaks close to each other
sels <- ohun::merge_overlaps(sels, pb = FALSE)
sels$dur <- sels$end - sels$start

tpks <- cbind(sels$start + (sels$dur - 6)/2, sels$peak)

# order by sample number
tpks <- tpks[order(tpks[, 1]), ]

# remove first peak tpks <- tpks[-1, ]

Then get the difference (in samples) between peaks (which is inversely related to the frequency)

sample_diffs <- diff(tpks[, 1])

Now get the start of each period in the sound clip

# # get envelope peaks amppk <-
# fpeaks(cbind(1:length(bout_1@left), abs(bout_1@left) /
# max(abs(bout_1@left))), plot = FALSE) # filter highest peaks
# amppk <- amppk[amppk[,2] > 0.95, ] # get the highest one
# first_max <- amppk[which.max(amppk[,2]),1] / bout_1@samp.rate

# plot
oscillo(bout_1, from = 0, to = 0.03, colwave = wave_col)

abline(v = cumsum(sample_diffs)/bout_1@samp.rate, col = "red", lty = 2,
    lwd = 2)

# plot
oscillo(bout_1, from = 0, to = 0.045, colwave = wave_col)

abline(v = cumsum(sample_diffs)/bout_1@samp.rate, col = "red", lty = 2,
    lwd = 2)

 

Takeaways

  • Amplitude cross-correlation seems to work fine for getting the periodicity although the position might be off a bit
  • Dynamic time warping does also a good job but is much slower
  • Dynamic time warping potentially more useful for shroeders in which the periods have variable lengths

 

 

Session information

## R version 4.1.0 (2021-05-18)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
## LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3
## 
## locale:
##  [1] LC_CTYPE=pt_BR.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=es_CR.UTF-8        LC_COLLATE=pt_BR.UTF-8    
##  [5] LC_MONETARY=es_CR.UTF-8    LC_MESSAGES=pt_BR.UTF-8   
##  [7] LC_PAPER=es_CR.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] viridis_0.6.2       viridisLite_0.4.1   warbleR_1.1.28     
##  [4] NatureSounds_1.0.4  seewave_2.2.0       tuneR_1.4.1        
##  [7] xaringanExtra_0.7.0 rprojroot_2.0.3     formatR_1.11       
## [10] knitr_1.42          kableExtra_1.3.4   
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.9        svglite_2.1.0     fftw_1.0-7        assertthat_0.2.1 
##  [5] digest_0.6.31     packrat_0.9.0     utf8_1.2.2        R6_2.5.1         
##  [9] Sim.DiffProc_4.8  signal_0.7-7      evaluate_0.20     highr_0.10       
## [13] httr_1.4.4        ggplot2_3.4.0     pillar_1.8.1      rlang_1.0.6      
## [17] uuid_1.1-0        rstudioapi_0.14   jquerylib_0.1.4   rmarkdown_2.20   
## [21] webshot_0.5.4     sketchy_1.0.2     stringr_1.5.0     igraph_1.3.5     
## [25] RCurl_1.98-1.9    munsell_0.5.0     proxy_0.4-27      Deriv_4.1.3      
## [29] compiler_4.1.0    xfun_0.36         pkgconfig_2.0.3   systemfonts_1.0.4
## [33] htmltools_0.5.4   tidyselect_1.2.0  gridExtra_2.3     tibble_3.1.8     
## [37] dtw_1.23-1        fansi_1.0.3       crayon_1.5.2      dplyr_1.0.10     
## [41] MASS_7.3-54       bitops_1.0-7      grid_4.1.0        DBI_1.1.1        
## [45] jsonlite_1.8.4    gtable_0.3.1      lifecycle_1.0.3   magrittr_2.0.3   
## [49] scales_1.2.1      cli_3.6.0         stringi_1.7.12    cachem_1.0.6     
## [53] pbapply_1.6-0     remotes_2.4.2     xml2_1.3.3        bslib_0.4.2      
## [57] vctrs_0.5.2       generics_0.1.3    ohun_0.1.0        rjson_0.2.21     
## [61] tools_4.1.0       glue_1.6.2        parallel_4.1.0    fastmap_1.1.0    
## [65] yaml_2.3.7        colorspace_2.0-3  rvest_1.0.3       sass_0.4.5